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SUMMARY 


This report presents the formulation, validation, and user’s manual for the computer 
program BL2D. The program is a fourth-order-accurate solution scheme for solving two- 
dimensional or axisymmetric boundary layers in speed regimes that range from low subsonic 
to hypersonic Mach numbers. A basic implementation of the transition zone and turbulence 
modeling is also included. The code is a result of many improvements made to the program 
VGBLP, which is described in NASA TM-83207 (February 1982), and can effectively 
supersede it. The code BL2D is designed to be modular, user-friendly, and portable to 
any machine with a standard fortran77 compiler. 

The report contains the new formulation adopted and the details of its implementation. 
Five validation cases are presented. A detailed user’s manual with the input format 
description and instructions for running the code is included. Adequate information is 
presented in the report to enable the user to modify or customize the code for specific 
applications. 
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«1> a 2, a 3 
a i,j 
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1.171 

hi 
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I 

i 

j 

k 

ki 

ke 
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L 

l 

h 

h 

M 
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N Pr 
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p 

p 

P 10 

Q 
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R 

Re 

r 

rf 

S 

s* 

T 


coefficients in the ^-direction differencing formula (eq. (18)) 
coefficients of the block-tridiagonal system given by eq. (49) 
diagonal elements of the linear system given by eq. (22) 
coefficients of the block-tridiagonal system given by eq. (49) 
superdiagonal elements of the linear system given by eq. (22) 
coefficients used in power law for viscosity (eq. (61)) 
skin-friction coefficient based on edge conditions, T w /{^p t u 2 e ) 

A0t/2 (eqs. (24) and (28)) 

ACf/12 (eqs. (24) and (28)) 
specific heat at constant pressure 

coefficients used in the block-tridiagonal system (eqs. (35), (41) and (47)) 
u I Ue 

an arbitrary function 

T/T e 

heat transfer coefficient, q* w /(T w - T aw ) 

H C . 

mdex m the stream wise (f ) direction 

flag; = 0 indicates two-dimensional flow; = 1 indicates axisymmetric flow 
index in the boundary-layer normal (£) direction 

value of index k at the edge of the inner part of the normal grid distribution 
index k that corresponds to the boundary-layer edge 
factor for stretching the grid in the C direction (eq. (62)) 

(pp/pePe) 
t 2l ll (eq. (1)) 
t 2 Hl (eq. (2)) 

Mach number 
sin 2 9 S 

Prandtl number, laminar value 

Prandtl number, turbulent value 

Stanton number based on edge conditions, h!{c v p\u \ ) 

pressure, dimensional 

pressure, nondimensional (eq. (53)) 

total pressure at boundary-layer edge (eq. (56)), nondimensional 

vector (eq. (14)) 

wall heat flux, dimensional 

gas constant, dimensional 

Reynolds number 

body radius 

residual, right-hand side of eq. (22) 
solution vector 

streamwise length along body surface, dimensional 
temperature, dimensional 
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t 

t 

t* 

U oo 
U 
W 
X 
z 
a 

P 

r 
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AC 

6 

c 

Os 

P 

t 

P 

Subscripts: 


aw 

e 

k = 1 
k = ke 
ref 
w 
0 

c.c 

oc 


transverse curvature term (used in eq. (1) and (2) only) 

temperature, nondimensional (eq. (55)) 

boundary-layer thickness, dimensional 

free-stream velocity, dimensional 

velocity in streamwise direction, nondimensional 

transformed normal velocity 

Cartesian coordinate in the body axis direction 

coordinate in the body-normal direction 

(i - ml 


streamwise intermittency distribution 
ratio of specific heats 
step size in C (eq- (15)) 
displacement thickness 
transformed normal coordinate 
shock-wave angle 
absolute viscosity 

transformed boundary-layer surface coordinate 
density 


adiabatic wall 
boundary-layer edge 
at the wall 

at boundary-layer edge 

reference value 

wall quantity 

stagnation condition 

partial derivatives in £ or C direction 

free-stream quantity 


Superscripts: 

* dimensional quantity 

' partial derivative in ( (except for x! , /, z*) 

n iteration number 

T transpose 

Abbreviations: 

BL2D boundary-layer two-dimensional program described in this report 

SI international (metric) system of units 

US U.S. customary or British system of units 

VGBLP boundary-layer program described in NASA TM-83207 

s/r subroutine 



1. INTRODUCTION 


This report presents the theory, formulation, and implementation of a fourth-order- 
accurate solution scheme for two-dimensional or axisymmetric boundary layers in speed 
regimes that range from low subsonic to hypersonic Mach numbers. A basic implementation 
of the transition zone and turbulence modeling is also included. The computer program that 
results from this work is available under the name BL2D. The basic formulation is based on 
the program VGBLP, which is described in NASA TM-83207 (ref. 1). 

The code BL2D is the result of many improvements made to the VGBLP program. The 
changes in formulation and implementation have been substantial; as a result, the coding 
for BL2D has been done from scratch. Thorough validation of the code has been completed; 
the author is confident that the code BL2D can effectively supersede VGBLP. 

The original theory and solution scheme presented in NASA TM-83207 is based on 
a second-order implicit scheme in the boundary layer normal direction. In contrast, the 
solution scheme in BL2D is based on a fourth-order-accurate compact Pade formula in the 
normal direction. This results in improved quality of the mean-flow profiles, which is a 
consideration of great importance for application to boundary-layer stability calculations. 
The Pade differencing scheme is discussed in detail in NASA CR-4531 (ref. 2) for three- 
dimensional boundary layers. 

The development of BL2D was also motivated by other considerations. The VGBLP 
program was written for machines of limited core memory and used many nonportable 
constructs, such as name lists and equivalence statements. It was also written in an 
older version of fortran and was designed for optimum performance on slower machines. 
Modification for specific applications was cumbersome, and updating the program with new 
models for turbulence and transition was difficult. 
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The present code is written in standard fortran77 and is portable (i.e., no machine- 
dependent features). The program package takes advantage of the user-friendly features of 
the UNIX operating system such as ‘make’ and file directories. Individual subroutines are 
now in separate files with ‘include’ statements that supply the common blocks. The code 
has been completely reprogrammed with emphasis on easy to understand and structured 
programming at the expense of some additional computation. An easy to use input format 
has been adopted. Output of solution quantities has been improved and is designed to be 
used for easy graphical interpretation. The code is easily modifiable to allow replacement of 
the built-in transition region and turbulence models with user-supplied routines in a modular 
form. The accuracy of the present method is independent of the normal grid distribution. 
With the computational speeds that are achievable today, the code can be run interactively 
on a desktop machine. 

The BL2D code has been validated with the original VGBLP test cases. The results are 
summarized in this report. A user’s manual is also included in this report. Sufficient details 
about the program have been provided in the manual so that it can be easily customized by 
the user for specific applications. The program package also includes the original VGBLP 
program files and input files for the five test cases described in reference 1. The complete 
program package and documentation are archived in the NASA Langley computer system 
and can be made available per individual request. 

The author would like to request that users communicate to him any errors, omissions, 
or desirable modifications to the report or the computer program. An updated description of 
such revisions will be maintained by the author and supplied with the software and report. 
The e-mail address of the author is v.iyer@larc.nasa.gov. 
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2. FORMULATION 


2.1 Equations in Vector Form 

We start with the transformed equations given in NASA TM-83207 (ref. 1). These 
are equations numbered (27) through (29) (page 15 of the report). A fourth-order-accurate 
Pade differencing scheme is now implemented, similar to the procedure outlined in NASA 
CR— 4531 (ref. 2). Let us use the following notations in place of the original notations used in 
NASA TM-83207 for the sake of simplicity, as well as for partial conformity with notations 
in NASA CR-4531. (See fig. 1.) 



— ►x* 

Figure 1. A sketch of notations used in boundary-layer formulation. 


We use w in place of V: H in place of 0] £ in place of tj for the transformed normal 
coordinate; and the prime symbol ' for the normal derivative (d/d(). The subscript £ is 
used to indicate the partial derivative in the £ direction. In addition, we use the following 
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notations: 


h = t 2j le ( 1 ) 

h = t 2j lt ( 2 ) 

L = F' (3) 

/ = H' (4) 

The equation set can then be rewritten as 

w — — F — 2 !fF c 

wF' - ( hF ')' = - 2 £FFs - p(F 2 - H ) ( 5 ) 

wH' - ( l 2 H ')' = - 2 (FHt + ah(F') 2 


Also note that in the equations above j = 0 or 1, depending on whether the flow is 
two-dimensional planar or axisymmetric; / = ( pfj,)/(p e fi e ); 1 is the eddy viscosity function 
^1 + l is the eddy viscosity function ; t is the transverse curvature 

term r/r w ; /? = and a = (7 - 1 )M 2 . 

To achieve a vectorial representation of the system prior to application of the fourth- 
order Pade differencing scheme, let us denote 


Q = 


( F \ 

Fw — l\L 

H 

\Hw - l 2 I ) 


( 6 ) 


The vector Q' can now be obtained from the system of equation set (5) as given below 
(with substitution for w'): 



L \ 

—F 2 (l +13)- 4 ZFFs + pH 

I 

—2^(FH)^ — FH + al\L 2 / 


( 7 ) 


4 


The element Q ^ of the equation above is obtained as 


Q' 2 = ( Fw - hL)' = Fw + wF' - (hF 1 )' = Fw - 2£FF^ -p (F 2 - H) 
= F{-F - 2 (F(.) - 2 tFFf: - P(F 2 - H) = -F 2 (l + P) - 4£FF* + pH 
The element Q'± of equation (7) is obtained as 


Q\ = ( Hw - l 2 I)' = Hw' + wH' - ( l 2 H ')' = Hw - 2 {FH^ + ah(F') 2 
= H (— F — 2 {Ft) - 2£F + ah (F') 2 = -2 t(FH\ — FH + ah (F') 2 
The second derivative Q n required in the differencing scheme is 

( L ' \ 

Q n = -2FL(l+P)-4Z(FLt + LFz)+pI 

\-2((FI + HL ); (FI + HL) + 2 a(hLL' + L 2 l[) ) 


( 9 ) 


(10) 


The variables L' and I' are obtained by differentiating the second and fourth elements of 
equation (6) with respect to ( and substituting for w' . The resulting expressions are 


l\V = F(2£Fjc + pF) + L(w- l[) - pH 


( 11 ) 


hi' = 2£FH{ - ah(F') 2 + l(w - l' 2 ) (12) 

The final form of Q" is, then, 

/ £ [F(2{Ft + pF) + L(w -l\)~ PH] 

q» = -2FL(\ +p)~ 4 £(FL f + LF{) + pi 

it 2tFHt-ahL 2 + I(w-l' 2 )\ 

\—2£(FI + HL )s - (FI + HL) + 2aL[F(2 ^ + PF) + Lw - pH] 

2.2 Discretization 

In accordance with the method in NASA CR-4531, we apply the fourth-order Pade 
differencing formula in the normal direction. This two-point compact scheme is defined in 
terms of the variable and its two higher derivatives. In the present case, if we assume that i 
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(i.e., the index in the surface direction £) remains constant and that k is the normal direction 
index, then the discretization at the midpoint of k and (k — 1) is written as 


Qt - Qt-i - % (<& + <&-.) + (<# - «*-•) + 0 ( A < 5 ) = 0 < 14 > 


ac = a - a-i 


(15) 


The differencing in the surface direction £ is accomplished to second order (or first order in 
some regions). For example, in the C direction, if we assume that the index k is fixed, then 

Qt = aiQi + {aQ} 

{aQ} = G2Qi-l + a iQi—2 


(16) 


The values of the coefficients a, are dependent on the location of the point in the streamwise 
direction. In the general case, the £ differencing is second-order accurate. In accordance 
with the parabolic nature of the equations, the £ derivative is obtained by the three-point 
upwind-differenced formula: 


(f()i — a ifi + °2/i-l + a zfi—2 

a, = (A£? - A5?.,)/A; a 2 = -Atf/A; a 3 = A&,/A (17) 

A = A£,A£i_i(A£i + A£^_i); A£ t — £» — £* — 1; A£j_i = £j-i — £t- 2 


When i = 2, the first-order formula with just two points is used, which results in the 
coefficients 

ai = 1/A£2 = l/(£2 - £1); 0,2 = — ai; <13 = 0 (18) 


A linear blending of the first- and second-order formulas can also be used in a specified 
region. (See section 3.3 “Streamwise Gradients” for details.) Thus, the short notation {aQ} 
used in equation (16) can be expanded in terms of the coefficients given above. 
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Jl 

We define a solution vector {5} = {F, L, H, 1} . Because the equations are nonlinear in 
{5}, Newton linearization is used to convert the system to a linear matrix inversion problem. 
If the superscript n denotes the current iteration stage, then we define {£5} as 

{<55} = 5” - 5 n_1 (19) 

A linear system is now set up to solve for {<55} in terms of the solution at iteration level 

r\ 

n — 1. For example, a term that involves (F") is written as 

(F n ) 2 = (F n_1 -f- <5F) 2 w (F n-1 ) 2 + 2F” -1 <5F (20) 

In the following equation set, the superscript n — 1 is dropped and is taken to imply the 
known values of {5} at iteration n — 1. A few examples of the linearized formulas with this 
notation are given below: 

(F n ) 2 = F 2 + 2F <5F 

F£ = Ft + ai 8F 

( 21 ) 

(Ftf)” = (FH)^ + ai (H6F + F6H) 

F n F£ = FF$ + 8F(aiF + F^) 

2.3 Linearized System 

The system is explicit in ( because of the choice of the finite-differencing scheme and is 
implicit in the surface- normal direction. The linearized system is represented at location i, 
which corresponds to the solution at iteration level n as 

’••• «{.*{«•••] {«?} = {r,‘J (22) 

where af m and bf m are elements of the (4 x 4) blocks in the diagonal and superdiagonal 
locations of the linearized block-bidiagonal system. The superscript k denotes that the 
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discretization corresponds to the midpoint of k and k — 1 points. The index / varies from 1 
to 4, depending on which element of equation (14) is being discretized; m varies from 1 to 
4, depending on which element of {65} it multiplies. The (4 x 4) blocks and m 

are the only nonzero blocks in the system above because of the two-point compact scheme 
in the k direction. The (4 x 1) vector {rf } corresponds to the residual of equation (14), 
based on the solution {5} at iteration n — 1. 

Let us write down the discretization of the first element of the system represented by 
equation (14), which can be written as 


(F) k ~ (F ) k _ i 


+ if {j;[ F ( 2 f F < + P F ) + L ( w ~ 'i) " W}, 
- ^{i[F(26F f + W) + L ( w - 4) - } t i = 0 


where 


ACj k = (k - Ct-1 


(24) 


We can now construct the elements of the blocks a* m j , > an< ^ { r f } f° r ^e case 

l — 1 from above by using the linearization procedure explained earlier. For example, 
the coefficient is the coefficient of (from the first element of eq. (11)), which is 

discretized at (k — ^), and is the corresponding coefficient of 6F k . 

Oil = “I ~ c 2Jfc( e ll)jfc-l 

ai2 = -ci* - C2*(ei2) Jt _ 1 

(25) 

«13 = — c 2Jb( e 13)fc_i 

ai4 = 0 
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6n = 1 + c 2 jt(en) fc 


&12 = ~Cih + c 2 *(e i 2 )jt 


^13 = c 2 *(ei3)* 


( 26 ) 


614 = 0 


r i — ~Fk + Fk-i + cik{Lk + Lk- 1) — C2k(qu)k + c 2 jt(<?n ) fc _ 1 


c\k 

C2k 



e n = l-[2f3F + 2t(a 1 F + Fz)] 
ei 2 = (w - l'i)/h 


e i3 = -0/h 


gn = Y [F(2£F( + f3F)+L(w- l[) - 0H] 


(27) 


(28) 


(29) 


(30) 


Similar expressions can be derived from the remaining three elements of the system that 


are represented by equations (6) — (8). The expressions for these elements of 


l l,m 


•i 




and {r*} are given below, preceded by the corresponding discretized equations. 


(Fw - l,L) t - (Fw - - ^-{-F 2 (l + 0) - 4 (FF ( + 0H} k 

~ {-f 2 (l + 0)- 4{FF ; + 0H} k _ k 

+ + « - i£{FLf + LF ( ) + 01} t 

~ ^{-2FI(1 +0)- 4 {(FL( + LF ( ) + 01}^, = 0 


(31) 
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a 2 i = ~Wk- 1 - CiJb(C2l)jb-l “ C2fc(e22)jt-l 


a 22 = (^l)jt-l “ c 2k( e 21 )fc— 1 

023 = ~ c lk0 

024 = ~C2k(3 

hi = w k - cut(e 2 i) fc + c 2 fc(e 22 )j t 

&22 = ( — ^l)jt + c 2*( e 2l)/fc 

hz = ~Clk0 
&24 = C 2 kfl 


( 32 ) 


(33) 


r 2 = — (?2l)jt + (92l)fc-l + Cl/t(922)jk + Clfc(^ 2 2)jfc_ 1 - C 2 Jb(?23)fc + C 2 fc(?23)*_i (34) 


e 2 i = — 2F(1 + 0) — + F{) 

e 22 = — 2L(1 + /?) — 4£(ai£, + L 


(35) 


g 2 i = Fw - l\L 

q 22 = -F 2 ( 1 + /?)- 4^FFjt + 0H (36) 

923 = -2FL(1 + 0)- 4Z{FLs + IF«) + 01 


fflrW-r + + 

- ^{^[ 2 « F77 < - “W 2 + 7 (" " ( 2 )l} t t = 0 



k 


(37) 
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031 = — C2Jfc(e3l)fc_i 

032 = c 2 it ( e 32 ) jb — 1 

033 = — 1 — C 2 fc(e 33)jt_i 

034 = ~Cn ~ C2Jt(e34)i_i 


( 38 ) 


^31 = C2fc(e3l)fc 

f>32 = C2k{eZ2)k 

633 = 1 + C 2 k(ezz) k 

634 = -Cik + C 2 jt(e34) fc 


(39) 


r 3 — ~Hk + Hk-i + cik(h + h-i) — C2*(93i)jfc + c 2 jt(93i)jt-i (40) 


«a = £[2 1«<] 

632 = — [— 2or/i£] 

‘2 

633 = f [2^0!] 

h 

ej4 = ( w - l' 2 )/l2 


(41) 


931 = 4 [2( /■'//, - aliL 2 + I(w - (')] 


(42) 
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( Hw - l 2 I) k - ( Hw - hI)k-\ 

- ^{-2 £(FH\ -FH + ahL 2 } k - ^{-2 ((FH\ - FH + *1,1?)^ 

+ ^*{—2 i(FI + HL) ( - (FI + HL) + 2aL[F(2^ + f3F) + Lw-/3H]} k (43) 

_ ^-{-2 '£(FI + HL\ - (FI + HL) + 2aL[F(2^ + 0F) + Lw - &H } } fci 

= 0 


041 = — c lJb( e 4l)ifc_l - c 2Jfc( e 42)jfc_i 

042 = -cu(e 43 ) fc _ 1 - C2fc(e 4 4) Jfc _ 1 

0 4 3 = -Wi fc-i - cu(e 4 5)jt_i - c 2 jfc(e 4 6)jt_i 

0 44 = - c 2 jfe(e 45 ) Jfc _ 1 


(44) 


641 = -ClJt(e 4 i) fc + C2k(e42)k 

6 42 = -c^e^)* +C 2 *(e 44 ) t 

643 = w k - cu^s)* + c^e^)* 

644 = (~h)k + ^(e^)* 


(45) 


J*4 = -(94l)fc + (94l)*_i + Cifc{942}jt + Clit{942}jk_i - C 2 jt{943}fc + ^2Jk {^43 } fc _i (46) 
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c 4 i = -2t(Ht + H ai )-H 

e 42 = -2£(7 e + /ai) - 7 + 2aL[2i{a x F + F e ) + 2/3 F] 
e43 = 2otl\L 

e 44 = —2 + ifai) -// + \awL + 2a[F(2£F,r + 0F) - /?//] 

645 = — 2^(Fjc + Fa\) — F 

646 = — 2^(F^ + -F01 ) — .7 — 2ot(3L 


( 47 ) 


941 = tfw — hi 


942 = -2£(F70* - FH + a/xl 2 ( 48 ) 

943 = — 2£(F7 + 7/7)^ - (F7 + 7/1) + 2a7[F(2^ + /9F) + 7u; - /?//] 


The resulting system is block tridiagonal of the order (4 x Jte x Are), where ke is the total 
number of points in the normal direction. The system is block tridiagonal because of the 
shift of rows that results from the implementation of the boundary conditions, as explained 
in references 2 and 3. Each row of the system at level k can be written as shown: 


f a l 1 

a 12 

a k 

a n 

a 14 \ 


*11 

*12 

*{ 3 

a 21 

a k 

a 22 

a k 

a 23 

a 24 


*21 

b 22 

*23 

0 

0 

0 

0 


4t l 

Ujf 1 

a^t 1 

a^ 1 

^0 

0 

0 

0 j 


/+ 1 
a 42 

a ii 1 

a 43 


*14 \ / 0 
*24 

„*+l 

a 34 

a *+i / 

a 44 / 


\b: 


0 

0 


0 

0 


0 \ 

0 


b^t 1 4 +1 4+ 1 b^ 1 

f+i b l 2 +i b l\ 1 b l\i 

°41 °42 °43 °44 


I 


6L k 

6H h 

{ bh ) 


> = < 



( 49 ) 
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The boundary conditions imposed at the wall and the boundary-layer edge modify the 
system. For the adiabatic wall condition, the first row of the block tridiagonal system is 
modified as 



( 1 

0 

0 

0 V 


0 

0 

0 

0 V 



l SFl ) 

f° ) 


0 

0 

0 

1 


0 

0 

0 

0 



* L i [ . 

0 1 


a !i 

a 32 

2 

a 33 

°34 


^31 

^32 

^33 

^34 


< 

8 Hi f 

r 3 | 


Wi 

a 42 

a 43 

“44/ 


V&41 

^42 

^43 

Ha) 



[8h ) 

UjJ 


Similarly, at the boundary-layer edge, the edge velocity and temperature are prescribed; 
this modifies the last row of the system as shown below: 


- 

(4’ 

ke 

a 12 

„ke 

a 13 




Lke 

°12 





' 8F ke ] 


ff] 


a 21 

a ke 

a 22 

n ke 

“23 

«24 



hke 

°22 

^23 




8L ) be 

> =2 4 

r k 2 e 


0 

0 

0 

0 


i 

0 

0 

0 



8H ke 


0 


Vo 

0 

0 

0 ) 



0 

1 

0 ) 



4 8Ike j 


0 J 


2.4 Inversion of the System 

At an y location i, assume that we have an initial estimate of the solution profiles F and 
H as a function of ( and their derivatives F' and H' . Also assume that we have an estimate 
of the transformed values of w and the viscosity ratios and assume that the edge 

conditions are known. With the knowledge of the upstream profiles, the system described 
by equation (49) is defined. Inversion of this block- tridiagonal system by standard routines 
yields the solution SF, 8L, SH and 81. However, this solution is based on the linearized 
formula. The solution is updated, and the calculation is repeated until convergence is 
achieved. 

Subsequent to convergence at a given station, the solution is advanced to the next step 
in the £ direction. The solution profiles at two upstream stations are saved at any given 
streamwise station. Lack of convergence is indicative of incipient boundary- layer separation. 
This can be confirmed by looking at the upstream history of C f e . 
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3. PROGRAM DESCRIPTION 


3.1 Edge Conditions from Inviscid Inputs 


The program can be run in SI or US units. The basic free-stream inputs are Moo, To, 
and To or , and Too- The isentropic formula 



( 52 ) 


supplies the missing quantities. The normalized pressure quantities are 


P 0 


Pq Poo 

Poo VI' Pc ° Poo Ul 


( 53 ) 


Note that pressure quantities are normalized by pooU^ and that normalized quantities are 
indicated by lower case. The gas properties Np r and 7 and the gas constant R are input. 
The reference temperature is defined as 


T te { = 


U 2 

^ OO 

h/h - 1)1^ 


( 54 ) 


The normalized total pressure and temperature values are, then, 

,,= 2 + F« ; = (7 - 1)1 Ml (55) 

To calculate boundary-layer edge velocity, the conditions on body surface must be 
calculated. If no shock is present between the free stream and the body, the normalized 
total pressure at the boundary-layer edge p 10 is equal to p 0 - If a shock is present, the 
condition behind the shock is calculated from the oblique shock-wave formula, based on the 
input shock-wave angle. If we assume for now that no shock curvature is present (hence, no 
variable entropy effects are present), then the computation requires the value of the shock 
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wave-angle 0 S as an input quantity. With Mq defined as sin 2 0 s , the normalized total 
pressure behind the shock is evaluated from 

f (7 + l)M» (7+1) (56) 

+ — (-r — l) J 


PlO 


The edge velocity normalized by the free-stream velocity is obtained as 


U e = 


\ 


2 4 _ © v } 


(57) 


where p e is the input normalized boundary-layer edge pressure. The normalized edge 
temperature is then 


t e = t o - 2«e 


(58) 


The edge density becomes 


pe — 


IPt 


(7 - l)*e 


(59) 


The edge viscosity p e (normalized by the reference viscosity p* ef evaluated at T re f) is 
calculated as a function of t e by using either the Sutherland law 

- ilL - #1-5 ( l±k.\ 

^ Pref 6 Ve + <r/ 


.5 


j»i 

p* et = 1.458 x 1° ~ & jrf+¥ Pa ’ SCC ° r 


jU.5 

! , = 2.27 x 1CT 8 ref 


Piet 


( lb • sec^ 


Tref + T r 

t r = Tr/T ld ] T r = 110.33 K or 198.6°R 


or the power law 


Pe = <? 


(60) 


pU f = 

Cl = 5.0231 x 10 -7 (SI units); ci = 7.1738 x 10 -9 (US units); c 2 = 0.647 


(61) 
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The program input is in the form of dimensional edge-pressure distribution at surface 
locations not necessarily coincident with the actual computation locations. The input 
quantities are interpolated to the required locations with a piecewise polynomial fit. 

3.2 Normal Grid Distribution 

Because the Pade formula is a compact scheme that is based on the solution variables 
and their derivatives at the two points that span the local cell center, a stretched grid can 
be employed without degradation of the fourth-order accuracy of the method. For laminar 
flows, no stretching is required in the transformed plane. However, for turbulent flows the 
large velocity gradients in the inner layer and relatively smaller gradients in the outer portion 
of the boundary layer require that some form of boundary-layer stretching be employed. 

The original code employed an exponential stretching in the normal direction such that 
the ratio k s is a constant. 

Ct+l ~ Ob _ , 

AC* ~ C* - Cfc-i “ s 

The normal grid distribution is then obtained as 

Cfc = Cke (^fcke - 1 (fc s > 1) (63) 

The input options to generate the normal grid in VGBLP were: specify k a , C* e , and ke 
(IGEOM = 1 in the program) or specify A (2 = C 2 — Ci> C*e, and ke (IGEOM = 2 in the 
program). This type of exponential distribution sometimes results in inadequate resolution 
near the edge of turbulent boundary layers because of the large C values required. This 
inadequate resolution in turn produces large truncation errors, especially when points are 
added in the normal direction to account for boundary-layer growth. 
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In the program BL2D, a two-part mesh distribution is implemented. (See fig. 2.) The 
mesh distribution in the inner part (delineated by 0 < C ^ C»> 1 < k < ^0 foN° ws 
the same exponential distribution as in VGBLP (eq. (63)). The input parameters are 
(i, k s , and ki (computer variable names ZI, AK, and NZl). The mesh distribution in the 
outer part is uniform in C, based on ( e , Ci, and ( ke - ki) (computer variable names 
ZMAX, ZI, and (NZ-NZI)). Obviously, the five input quantities ( k s , Ci, ki, (e, and ke) 
must be carefully chosen to avoid abrupt changes in the mesh step size at k ki. The 
program gridcheck . f can be used to arrive at the optimum distribution before running 
the boundary-layer program. To set a single exponential stretching for the entire region, set 
Ci = Ce and ke = ki. The advantage of this two-part distribution is evident when points 
must be added to accommodate boundary-layer growth. 


10 

8 

6 I 

4 

2 


C = Ce k = ke 


Outer distribution 


C_=k 


Inner distribution o0 °' 


C e = 8.0 , C/ = 4.0 
ke = 45 , ki = 30 , k s = 1.05 

, I L , . df I 


0.0 


0.1 0.2 

Mesh spacing, 


0.3 


Figure 2. Two-part mesh distribution used in program. 

In the event that the user wants to implement a different normal grid distribution, this 
can be set up easily in subroutine grid. The addition of points to cover boundary-layer 
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growth is done in subroutine add. This addition of points is also easily modifiable if the 
user wants to change the present logic (i.e., each new point is added at a constant step size). 

3.3 Streamwise Gradients 

The streamwise grid distribution is specified in the input by means of the number of 
steps and the corresponding step sizes As*. The streamwise gradient is calculated based on 
a second-order formula that can be blended into a first-order formula for a certain region. 
The transformed streamwise coordinate £ is 

3 

( = J PeUeVe(r w ) 2j ds ( 64 ) 

0 

The £ derivative from a three-point upwind-differenced formula (valid for i > iord2) is 

(/(),• = <*ifi + 02/1-1 + 03/1-2 

«l = (A£, 2 - A&O/A; a 2 = -Af?/A; a 3 = A^Lj/A (65) 

A = A&A£,-i(A& + Afc-i ); A 6 = 6 - &_ i; A&_i = 6-i - 6-2 

When i < iordl, the first-order formula with just two points is used, which results 
in the coefficients 

a\ 1/A6 = 1/(6 - 6); °2 = — «i; «3 = o (66) 

For iordl < * < iord2, the first-order formula is blended to the second-order formula 
with a linear variation. 

3.4 Starting Solutions 

Depending on the geometry, the starting solution is generated from the similarity 
assumption (for a body with a sharp tip or leading edge) or from the equations that are 
valid for a two-dimensional stagnation point (for a body with a blunt tip or leading edge). 
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The similarity solution is obtained by dropping all streamwise gradient terms. In the 
program, this step is done by setting the coefficients a\, a. 2 , and <23 to 0. 

3.5 Implementation of Boundary Conditions 

Mass injection or suction at the wall is input into the code in the form of the normal 
momentum p* w Ww (in dimensional units, (Pa-s)/m or (lb-s)/ft 3 ). These values are first 
interpolated to the boundary-layer grid locations from the input distribution. The value 
of the transformed normal velocity w at the wall is obtained from the relation 


w 1 = 


p e rw(p* e u*) 


The continuity equation (eq. 5, line 1) is integrated with the trapezoidal rule; the initial 
value is specified as shown above. A negative value for Pw w w indicates suction. 

The energy boundary condition at the wall is specified in the input by selection of the 
variable iwall. A value of 0 indicates adiabatic wall conditions, and a value of 1 indicates 
that the wall temperature T w is to be specified. A value of 2 indicates that the dimensional 
heat flux at the wall q w * in dimensional units (W/m 2 or Btu/ft 2 -sec) is to be specified. 
In this case, the temperature gradient at the wall is specified in terms of the transformed 
variable Hj as 


(ulC t) 


— Qw — 


pePe^e^e^tv 

^TiV^Npr 


Note that the input and output of heat flux in US units is in Btu/ft 2 -sec. Within the 
program, the units of lb/(ft-sec) is used (1 Btu/ft 2 -sec = 778.26 lb/(ft*sec)). 

3.6 Variable Entropy Formulation 

In the case of highly curved shocks in front of blunt bodies in high-speed flow, the change 
in edge conditions that results from the entropy layer may need to be taken into account. 
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A detailed discussion of the approach followed here is given in references 4 and 5. Briefly, 
this step involves an initial complete boundary-layer calculation that neglects this effect. 
Subsequently, from the boundary-layer solution information a mass balance is performed 
at each streamwise station to compute the radius of the corresponding stream tube at the 
shock location. The local shock angle and the corresponding post-shock edge conditions are 
computed from this calculation. The boundary-layer calculations are then repeated with the 
new edge conditions until the variable entropy convergence criterion is satisfied. 

3.7 Transition Zone Modeling 

The location of transition is either explicitly specified or calculated from the stability 
index parameter reaching a specified value. For details, refer to reference 3. This simple 
representation of transition onset can be replaced by a user-supplied model. The extent of 
the transition zone is calculated in the program by using a simple correlation that is based 
on the local Reynolds number or it can be explicitly specified. See references 1 and 4 for 
a complete description. The intermittency function T is used to characterize the transition 
zone and is calculated as a function of £. AT value of 1 signifies the beginning of the fully 
turbulent zone. 

3.8 Turbulence Modeling 

The algebraic turbulence models used in BL2D are described in references 1 and 4. 
Two choices are available: a general mixing length model (KODVIS = 1); or a two-layer 
eddy-viscosity model (KODVIS = 2). The intermittency factor, T, is multiplied by the e//i 
values and used to model the transitional zone. The program is designed to allow a user to 
implement any desired model, for example, two-equation closure models. 

3.9 Internal Flows 

These flows are treated in the same manner as external flows. The curvature terms will, 
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however, involve a change in the sign. If curvature terms are neglected, then the solution is 
identical to the external flow. If the curvature terms become significant, this indicates that 
the boundary-layer thickness values are significant enough to alter the inviscid pressures for 
an internal flow. This in turn indicates significant viscous-inviscid interaction. In view of 
this, for internal flow cases, computation should be done neglecting the curvature terms. 

3.10 Output of Physical Quantities 

The output of wall quantities and boundary-layer integrated parameters are determined 
based on user-selected codes. Solution profiles at any station can be output by enabling the 
corresponding flags. Details are given in chapter 5, entitled “User’s Manual.” 
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4. VALIDATION 


Five cases are presented here to validate BL2D with VGBLP. These cases are identical 
to those presented in reference 1. Most capabilities of the code are covered between these 
cases. Because the results compare well, the author is confident that BL2D can effectively 
supersede VGBLP. 

4.1 Case 1 : Flow Past a Flat Plate at Mach 2.8 

A sketch and description of the flow case and conditions are given in figure 3 and table 
1. As indicated in the sketch, the flow has laminar, transitional, and turbulent regions. The 
thermal boundary condition at the wall is adiabatic. 

The input boundary-layer edge conditions are constant and are specified in terms of two 
data points. The step size is variable and the solution is obtained at 62 stations. 

The results from BL2D are compared with the VGBLP results for identical input 
conditions in figures 4—9. Figure 4 is a comparison of the skin-friction coefficient based 
on the edge conditions (C/.e), 99.5-percent boundary-layer thickness (f*), and displacement 
thickness (6*). Figure 5 is a close-up of the same results around the transition region 0 < 
x < 0.04 in. These parameters are dependent on the wall gradients and the shape of the 
profiles. Excellent agreement has been obtained. 

Figure 6 shows a comparison of the wall temperature ratio ( T w /To ) and the streamwise 
intermittency parameter. Figures 7—9 are comparisons of the solution profiles F, H, and w 
at six locations in the flow (four in the turbulent region and one each in the laminar and 
transitional regions). Again, excellent agreement has been demonstrated. 
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Figure 3. Case 1: geometry and conditions. 



Table 1. Case 1: Description and Input Summary 


Description 

Case 1 : Supersonic Flow over Flat Plate 
Free-stream conditions 

M-inf =2.8 

P-tot = 4.14 x 1CT6 N/m 2 
T-tot = 311 deg K 

Wall conditions 

Adiabatic wall (IWALL = 0) 

No mass injection 

Flow type 


No shock (WAVE = 0 deg) 

2D (J = 0) 

No stagnation point ( IBODY = 2) 
Constant entropy ( IENTRO = 1) 

Viscous terms 


Laminar 

Transitional (SMXTR = 2400); for s > 0.005 m 
Turbulent (TLNGTH = 2); for s > 0.01 m 

Other 


Body opening angle PHII = 0 

Solution for 0 < s < 0.25 m at 62 stations 
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Symbols refer to VGBLP solution 
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Figure 4. Case 1 : comparison of BL2D and VGBLP solutions (C fe , 
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Figure 7. Case 1 : comparison of BL2D and VGB 
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Figure 8. Case 1 : comparison of BL2D and VGBLP solution profiles (T/T e ). 
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Figure 9. Case 1 : comparison of B 


4.2 Case 2: Flow Past a Waisted Body at Mach 1.7 


This case involves a flow with an oblique shock that results from a sharp-tipped ax- 
isymmetric body; the resulting boundary layer has a pressure gradient of varying sign and 
magnitude. Figure 10 shows a sketch of the body and the variation of the surface pressure. 
The body has a short laminar region followed by a short transition region. The onset and 
the extent of transition are explicitly specified. Table 2 presents the flow and input details. 

Figure 11 shows a comparison of the skin-friction coefficient based on edge conditions 
(C/ >e ), 99.5-percent boundary-layer thickness (t*), and displacement thickness (6 ). Figure 
12 shows a comparison of the wall temperature ratio ( T w /To ) and the streamwise intermit- 
tency parameter. Figures 13-15 axe comparisons of the solution profiles F, H, and w at five 
locations in the flow (four in the turbulent region and one in the laminar region). Again, 
excellent agreement has been obtained. 

For this slender axisymmetric body with appreciable boundary-layer thickness values in 
comparison with the body radius, the curvature terms are likely to be significant. The results 
presented above include this effect by enabling the curvature terms in the equation (option 
IW = 1). The results of a run made with the curvature terms disabled (IW = 0) are shown in 
figure 16. Previous results with the curvature effects included are also shown for comparison. 
A noticeable difference is evident in the solution, especially at downstream locations at which 
the boundary- layer thickness is significant. The present results also compare well with the 
VGBLP results obtained with the curvature terms zeroed out. 
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Figure 10. Case 2: geometry and conditions. 




Table 2. Case 2: Description and Input Summary 


Description 

Case 2: Supersonic Flow over Waisted Afterbody 
Free-stream conditions 


M-inf =1.7 

P-tot = 47511 N/m2 

T-tot = 297.8 deg K 

Wall conditions 


Adiabatic wall { IWALL - 0) 
No mass injection 

Flow type 


Oblique shock (WAVE = 43.523 deg) 
Axisymmetric (J = 1) 

No stagnation point ( I BODY = 2) 
Constant entropy ( IENTRO = 1) 

Viscous terms 


Laminar 

Transitional (SST = 0.0901); for s > 0.09 m 
Turbulent (TLNGTH = 2); for s > 0.18 m 

Other 


Body opening angle PHI I = 20 

Solution for 0 < s < 1.53 m at 165 stations 

Curvature term on or off (IW = 1 or IW = 0) 
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Symbols refer to VGBLP solution 
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Figure 12. Case 2: comparison of BL2D and VGBLP solutions (T, (TJTJ). 
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Figure 13. Case 2: com] 
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Symbols refer to VGBLP solution 
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Figure 16. Case 2: comparison of BL2D and VGBLP solutions (C f t , t ,8 ) with 
effect of zeroing out curvature term ( IW = 0 ) . 



4.3 Case 3: Flow Past a Sharp Cone at Mach 7.4 


This case involves a sharp cone with a half angle of 5° and a cooled wall. This case 
demonstrates the use of the wall boundary conditions in the code. Three subcases correspond 
to the prescribed conditions of no suction, suction, and blowing (cases 3.1, 3.2, and 3.3, 
respectively). The boundary-layer grid is clustered near the region in which the wall mass 
transfer has a step change. In all three runs, the wall temperature is a constant value. The 
flow is laminar. Details of the geometry and the input conditions are given in figure 17 
and table 3. 

Case 3.1: No Suction . Figure 18 shows a comparison of the skin-friction coefficient based 
on edge conditions (C/ >e ), 99.5-percent boundary-layer thickness (t*), and displacement 
thickness (S*). Figure 19 shows a comparison of the wall heat transfer rate and the 
Stanton number based on the edge conditions. Figures 20—21 are comparisons of the solution 
profiles F and H at five locations in the flow. The BL2D profiles are slightly fuller than the 
VGBLP profiles, probably due to the different numerics of the differencing schemes. 

Case 3.2: Boundary- Layer Suction . A constant suction rate of —0.090117 Pa-s/m is 
applied for s* > 0.096012 m. The resulting variations of C/ >e , t*, and 6* are shown in figure 
22. The variation of and the Stanton number axe shown in figure 23. The agreement 
between VGBLP and BL2D is good. The profiles of u/u t and T/T e are shown in figures 
24—25. The profiles are fuller as a result of suction. 

Case 3.3: Boundary-Layer Blowing . A constant blowing rate of 0.090117 Pa-s/m is 
applied for s* > 0.096012 m. The boundary layer separates at around s* = 0.11 m. The 
resulting variations of Cf >e , t*, and 6* are shown in figure 26. The profiles of u/u € and T/T e 
at three locations are shown in figures 27-28. Because of its proximity to the separation 
point, profiles at the third location of BL2D and VGBLP are slightly different. 
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Figure 17. Case 3: geometry and conditions. 



Table 3. Case 3: Description and Input Summary 


Description 

Case 3: Hypersonic Flow over Sharp Cone 
Free-stream conditions 

M-inf = 7.4 

P-tot = 4.14 x 10^6 N/m 2 
T-tot = 833 deg K 

Wall conditions 


Wall temperature specified ( IWALL = 1) 
No mass injection (Case 3.1) 
Prescribed suction (Case 3.2) 
Prescribed blowing (Case 3.3) 

Flow type 


Oblique shock (WAVE = 9.214 deg) 
Axisymmetric (J = 1) 

No stagnation point ( I BODY = 2) 
Constant entropy (IENTRO = 1) 

Viscous terms 


Laminar 

Other 


Body opening angle PH I I = 5 

Solution for 0 < s < 0.3 m at 85 stations 
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and VGBLP solutions (C, 




Symbols refer to VGBLP solution 
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Figure 19. Case 3.1: comparison of BL2D and VGBLP solutions (q, 
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Figure 
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Figure 2 ‘ 
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Figure 
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: comparison of BL2D and VGBLP soli 


4.4 Case 4: Hypersonic Flow Past a Blunt Cone 

Here, we look at the flow past a blunt cone at a Mach number of 20.3 in Helium gas (in 
which 7, R y and Np r are different from those in air). Further, the shock curvature in this 
problem indicates the need to use the variable entropy option IENTRO = 2 in the code. 
Figure 29 and table 4 give a description of the flow and input conditions. 

Figure 30 presents the results from BL2D and VGBLP that correspond to the first pass 
in the variable entropy iteration (ITE = 1). This computation is equivalent to running the 
code with no variable entropy iteration (IENTRO = 1). Good agreement is obtained for 
Cf,ei t 5 and S . Note that in this case the BL2D computation uses first-order streamwise 
derivatives because the edge conditions are not smooth. 

Figure 31 shows the same results at the end of the variable entropy iteration convergence 
(ITE = 3 in this case for both codes). Note that the calculation now involves the shock 
curvature determined locally from the slope of the shock front. Because the original shock 
shape and gradient data from reference 1 are not smooth, a spline smoothing was done 
(external to the program), and the smoothed coordinates and slopes were used in the BL2D 
and VGBLP computations. Results in figure 31 show that the variable entropy has a 
significant effect on the results. The two codes agree fairly well in the variable entropy 
mode; the slight difference in results near the stagnation point was traced to an error in 
initialization for ITE > 1 in VGBLP. Figure 32 shows a comparison of the heating rate 
at the wall. The rough edge data and the nearly separated boundary layer for s* > 0.02 m 
cause the slight difference in the results. 
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Figure 29. Case 4: geometry and conditions. 




Table 4. Case 4: Description and Input Summary 


Description 


Case 4: Hypersonic Flow over Blunt Body 


Free-stream conditions 


M-inf =20.3 

P-tot = 7 x 10 A 6 N/m 2 

T-tot = 289 deg K 

Wall conditions 


Wall temperature specified (IWALL = 1) 
No mass injection 

Flow type 


Curved shock, normal at stagnation point (WAVE = 90 deg) 
Axisymmetric (J = 1) 

Stagnation point ( IBODY = 1) 

Variable entropy (IENTRO = 2) 

Viscous terms 


Laminar, power law for viscosity 
Other 


Body opening angle PHII = 90 

Solution for 0 < s < 0.03 m at 121 stations 

Helium gas (PRL = 0.688, GAM = 1.6667, RSTAR = 2079.0 m A 2/ (sec^ degK) 
Shock wave coordinates input for variable entropy calculation 
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Symbols refer to VGBLP solution 
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variable-entropy results for iteration 1 . 
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Symbols refer to VGBLP solution 
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Figure 32. Case 4: comparison of wall heat flux (with variable entropy iteration enabled). 



4.5 Case 5: Turbulent Flow In a Convergent-Divergent Nozzle 

The free-stream or inflow Mach number here is 0.012058; the flow is rapidly expanded 
in the nozzle to an exit plane Mach number of 5. The flow is turbulent almost from the 
beginning. The flow and geometry parameters are given in figure 33 and table 5. 

A comparison of the displacement- thickness values from BL2D and VGBLP is given in 
figure 34. A log scale is used for 8 because of its exponential growth in the divergent portion 
of the nozzle. Figure 35 shows the standard plot of C/ >e , t *, and 8* as a function of s*. 
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5: geometry and conditions. 



Table 5. Case 5: Description and Input Summary 


Description 

Case 5: Flow in a Convergent-Divergent Nozzle 

Free-stream conditions 

M-inf = 0.012058 
P-tot = 3.45 x 10 A 6 N/m2 
T-tot = 377 deg K 

Wall conditions 

Adiabatic wall ( IWALL « 0) 

No mass injection 

Flow type 
No shock 

Axisymmetric (J = 1) 

No stagnation point { IBODY = 2) 

Viscous terms 

Turbulent 

Other 


Body opening angle PHI I = 0 

Solution for 0 < s < 0.62 m at 101 stations 
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values. 






5. USER’S MANUAL 


5.1 BL2D Program Structure 

The program is located in several directories. Please refer to figure 36 for the locations 
of various files in the program package. 

The subdirectory bl2d/source contains the Fortran files that correspond to the 
subroutines of BL2D. Herein can also be found a Makefile to automatically compile the 
required subroutines and create the executable code. Details for running the program are 
given in section 5.3, “Running the Program BL2D”. 

The subdirectory bl2d/inputs contains input files for the five test cases used in 
validation. (See chapter 4, entitled “Validation”.) This subdirectory also contains versions 
of the include files used for the test cases. These files contain the common blocks used in 
the source routines. Array dimensions are set with parameter statements. 

The subdirectory bl2d/lib contains some useful programs for the input check. Please 
refer to the Readme file in this subdirectory for further details. 

The subdirectory bl2d/vgblp contains a version of VGBLP used in the validation 
runs, along with corresponding input files. Many modifications have been incorporated in 
the original program to run on UNIX machines with a fortran77 compiler. Note that the 
original program statements are in upper case and the modifications that are incorporated 
are in lower case. The Readme file in this subdirectory explains how to run the test cases. 

The subdirectory bl2d/doc contains documentation on BL2D. Future updates will also 
be documented in this area in an Update . info file. 

A synopsis of the call sequence of various routines in BL2D and the solution logic is 
given in appendix A. Appendix A also contains a brief description of the purpose of each 
subroutine. 
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/vgblp 
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Figure 36. Program files structure. 









5.2 Description of Input File for BL2D 

A sample input file for BL2D (for Case 1) is shown below. This input file is the same as 
bl2d/inputs/inp . 1 included in the BL2D program package. 


VALIDATION OF BL2D WITH 

VGBLP CASE 1 MAi 

IUNIT 


AMACH 


PTS or ! 

1 


2.8 


0 . 414e7 

GAM 


IGAS 


I WALL 

1.4 


1 


0 

I BODY 


WAVE 


PHII 

2 


0.0 


0.0 

ZMAX 


ZI 


NZI 

40.0 


40.00 


81 

DFPTOL 


DHPTOL 


I ACC 

l.e-03 


l.e-03 


4 

DFE 


DHE 


IADD 

0.0001 


0.0001 


0 

IYINT 


KODAMP 


KODPRT 

2 


2 


1 

SMXTR 


SST 


TLNGTH 

2400.0 


l.e8 


2.0 

RSTAR 


IORD1 


IORD2 

286.96 


2 


3 

STEP SIZES (NX1 

VALUES) 



.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.001 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.00 5 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.005 

.010 

.010 

.010 

.010 

.010 





WALL OR 

PROFILE 

PRINT FLAGS (0, 

1 or 2), 

1 

1 

1 

1 

2 

1 

1 

1 

2 

1 

1 

1 

1 

1 

1 

1 

1 

1 

2 

1 

1 

1 

1 

1 

2 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

2 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

2 

1 





OUTPUT 

FOR PLOTS, NUMBER OF VARIABLES 1 

8 

34 

2 

25 

46 

NUMBER 


L (INVISCID INPUTS) 

2 


1 



XE 

RADE 

SE 

PESE 

TWSE 

0.000 

1.0 

0.0 

152552. 

. 0.0 

1.000 

1.0 

0.100 

152552. 

. 0.0 


2.8 FLAT PLATE, X-0.0 to 0.25 


TTS or TFS 
311.0 
J 
0 

IENTRO 

1 

AK 

1.10 

ITMAX 

100 

VELEDG 
0.995 
ROD VIS 
2 

PRT 

0.95 

ITEMAX 

1 


NX1 VALUES 


IFS (- for FS) 
1 

1 FT 
1 

CONVE 

l.e-02 

N2 

81 

IW 

0 

NX1 

61 

KTCOD 

2 

NXLIM 

62 

PRL 

0.72 


QESE 

0.0 

0.0 


39 


WWSE 

0.0 

0.0 


45 


40 


Most of the input parameters correspond to VGBLP. (In most cases, the same variable 
names are used.) However, the input format has been improved. The namelist input 
format is no longer used. Line headers are used in the input file to facilitate input in 
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an unambiguous and error-free manner. Common input errors, such as improper array 
dimensions and incompatible values, are trapped in the program. 

The input files for running the five validation cases are available in the bl2d/inputs 
subdirectory. These input files may be useful to the user in the selection of appropriate 
input parameters for a new case. 

A detailed description of the input file format is given in appendix B. (The same 
information is also available in the file bl2d/doc/inp . doc.) An alphabetical list of BL2D 
input variables with a short description is given in appendix C. In this section, some of the 
new features in the input are described. 

The free-stream conditions can be specified in terms of stagnation or static values (i.e., 
(PTS, TTS) or (PFS, TFS) ). The choice is indicated by the input variable, IFS; a negative 
value indicates that static values are to be input. 

As given previously, the normal grid distribution is input such that an inner and an 
outer distribution can be specified; a single exponential distribution can also be specified, as 
in VGBLP. See section 3.2 entitled, “Normal Grid Distribution” for more details. Also note 
that the program gridcheck. f in the subdirectory bl2d/lib can be used to arrive at 
the optimum distribution before the boundary-layer program is run. 

In BL2D, the solution convergence is tested on the maximum values of both dF' and 
dH' . The corresponding tolerance values are input. If the convergence tolerances are not met 
in ITMAX iterations, then a warning of no adequate convergence is issued and the solution 
march is continued. 

The flag iadd is a new feature by which the normal grid extent and dimension can be 
increased to accommodate the growth of the boundary layer in transformed variable space, 
if needed. Normal grid points are presently added at equal increments. If no more points 
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can be added due to a limit on the array dimensions, an error message will be printed. 
The criterion for adding normal grid points is based on sum of the gradient F' at the 10 
outermost grid points. If this sum exceeds the input tolerance value DFE and if iadd = 
1, the program adds points automatically. A similar test is done on the temperature profile 
gradient, based on the input value of DHE. 

The input of nxlim < nx enables the solution march to be stopped at an earlier step 
for diagnostic runs. 

Wall parameters or solution profiles are printed via flags. A value of 0 indicates no print, 
1 indicates that the wall parameters will print, and 2 indicates that both the wall parameters 
and the profiles will print. This flag is specified at the end of each marching step. A flag 
of 1 at the 10th location indicates that wall parameters at i = 11 (i.e., end of step 10) will 
be printed. Wall parameters are output to fort . 2, and profiles are output to fort . 8. In 
addition, up to 49 wall parameters can be output to file fort . 7 at all i locations (i > 1) 
for plotting purposes. A number code which ranges from 1 to 49 is used for the variables 
that can be output. The list of these variables is given in appendix D. Up to 10 variables 
can be output at a time per run. 

5.3 Running the Program BL2D 

The BL2D program is run in the subdirectory bl2d/source. The program expects 
input under the file name inp . dat and an include file that contains common blocks 
under the file name com. The maximum dimensions of the array are set in the file com via 
parameter statements. The array dimensions are as follows: 


nzm = maximum value of normal grid points 
nxm = maximum value of streamwise grid points 
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nmax 


greater of (nzm,nxm) 


num = number of inviscid data points 

nsm = number of points that define the curved shock 

set to 1, if variable entropy option not used 

numblm = number of points that define PRTAR array, 1 if kodprt 
= 1 or 2 

Case 1 example: 

parameter (nzm = 81,nxm = 62, nmax = 81, num = 2, nsm = 1, numblm = 1) 


To run Case 1, for example, copy com . 1 and inp . 1 (from subdirectory bl2d/inputs) 
to com and inp . dat, respectively; remove all . o files (only if a new com file is being used, 
as in this example); type make to create the executable code; and type a . out to run the 
program. Detailed output is written to the file fort. 2. Other outputs are written to 
fort . 7 and fort . 8. The input and include files for other cases are available as inp . n 
and com.n in the subdirectory bl2d/ inputs (where n refers to the case number). 

5.4 Modifying the Code 

Inevitably, modifications must be added to the code for specific applications. The 
modular structure of the code enables the user to perform this easily. Most often, additional 
output statements are added or some new parameters are computed. New transition 
zone models and turbulence models can also be easily added. The modifications can be 
incorporated in main . f by additional call statements or write statements. To assist the 
user in possible modifications, a list of the main variables used in the program and their 
definitions are given in appendix E. 
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Finally, the author requests that he be kept informed of details of errors, omissions, 
and suggested modifications to the report or the computer program, if any. An updated 
description of such revisions will be maintained by the author and supplied with the 
software and report (file bl2d/doc/Update . inf o). The e-mail address of the author 
is, v. iyer01arc.nasa.gov. 
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APPENDIX A. BL2D PROGRAM LOGIC 



□ Read input data (free-stream parameters, solution options and parameters, inviscid 
edge parameters) from file inp . dat 

□ Inputs echo-printed to fort. 2 


Subroutine:; 


□ Based on input free-stream conditions, calculate reference conditions and normalizing 
quantities 

□ Summary printed to fort. 2 




Subroutine • ■ 

□ From input step sizes ss (i) , calculate surface arc length s (i) 

□ From input edge data, interpolate for edge quantities at s ( i ) locations 

□ Computes normalized values of ue, te, tw , amue , roe; calculate £ 

□ Interpolated edge values written to fort. 2 


Subroutine* .\grid 


□ Generate initial normal grid based on input selections 

□ Normal grid in £ direction written to fort. 2 

y 

_»» « Entry point (1) for iteration for flows with curved shocks (variable entropy option); 
I edge velocity u e calculated based on current solution; ite is the iteration counter 
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'Subroutine ini t i 

□ Initialize profiles of F,H, F' ,w 

□ Set to zero H'\ the upstream profiles of F,H,F',H' ; set e = 1, l = l/Np r 

□ Call subroutine updte to initialize profiles of fi, /, / 1? l 2 , 

ir 

Entry point (2) : start do loop for solution march in the i direction (i = 1, nxlim) 
| Subroutine edge 


□ Calculate edge coefficients for the current i 

□ Evaluate M e , a, p 

□ Set wall boundary condition: if iwall = 0, H[ = 0, if iwall = 1, H\ = T w /T e \ 
ifiwall = 2, update of H[ is done in subroutine updte 

□ Special considerations for i = 1 similarity flow or stagnation point flow 


-*• • Entry point (3) for calculation with new normal grid 

I 

Entry point (4) for iteration loop; it is the iteration counter 


1 


Subroutine blk 


□ Solve momentum and energy equations with current known profiles 

□ Update F, H , F', H' 
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□ Obtain current dimensional profiles 

□ Calculate current boundary- layer parameters (required in subroutines trans and 
turb) 


Subrouting brans f | 

□ Basic transition onset model and transition zone model 

□ r calculated as T (0,1) 


V 

Subroutine turb • • •• ' 

□ Basic algebraic turbulence model (called if T is > 0) 

□ e, e calculated 

□ In main program update viscous terms again (call updte) if s/r turb was called 

7 

; ' Subroutine ccmvVii: 

□ Test for convergence based on maximum values of dF',dH' 


< ■ ■« Branch to point (4) if convergence criteria are not met 


Subroutine a&d 


□ Check for adding points in the normal direction 
_ . « If yes, call s/r updte and branch to point (3) 


MAIN PROGRAM 

] L. Recalculate dimensional profiles (s/r phys) ; call s/r nextep (save previous profiles) 
^ and branch to point (2) . 

— • Branch to (1) if variable entropy convergence criterion is not met 


f 

llllliil: 
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APPENDIX B. INPUT FILE FORMAT FOR BL2D 


Line 1 
TITLE 

Line 2 
IUNIT 

Line 3 
Values 
IUNIT 

AMACH 

PTS 

PFS 

TTS 

TFS 

IFS 


Line 4 
GAM 

Line 5 
Values 
GAM 
IGAS 


I WALL 


J 

IFT 

Line 6 
I BODY 

Line 7 
Values 
IBODY 

WAVE 

PHII 

IENTRO 

CONVE 


(no more than 79 characters long) 


AMACH PTS or PFS TTS or TFS IFS 


corresponding to line headings above 

0 for British units 

1 for SI units 
mach number 

free-stream total pressure (lb/ft2 or N/m2) 

free-stream static pressure (lb/ft2 or N/m2) 

free-stream total temperature (deg R or deg K) 

free-stream static temperature (deg R or deg K) 

negative value indicates that input values are PFS and TFS 

positive value indicates that input values are PTS and TTS 


IGAS IWALL J 


IFT 


corresponding to line headings above 
ratio of specific heats 

1 for Sutherland law for viscosity 

2 for Power law 

default values for the corresponding coefficients 
(vislcl, vislc2 for igas=l; vis2cl / vis2c2 for igas=2) 
are specified in s/r ref 
wall boundary condition; 

0 for adiabatic 

1 for wall temperature specification 

2 for wall heat flux specification (if other than zero) 

0 for 2-D 

1 for axisymmetric 

0 for locally similar solution (for checking purposes) 

1 non-similar solution (recommended) 


WAVE PHII IENTRO CONVE 


corresponding to line headings above 

1 for flow with a stagnation point 

2 for flow without stagnation point 

shock wave angle at x=0 (input 0 if no shock) 
opening angle of body at x=0, deg. 

1 for constant entropy 

2 for variable entropy 

convergence criterion for variable entropy iteration 
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Line 8 
ZMAX 


ZI 


NZI 


AK 


NZ 


Line 9 

Values corresponding to line headings above (normal grid parameters) 

- . , /- a T 


ZMAX 
ZI 
NZI 
AK 
NZ 

Line 10 
DFPTOL 


the maximum value of the transformed normal co-ordinate 
normal co-ordinate value at the end of inner distribution 
number of mesh points in the inner distribution 

stretching parameter for inner distribution (AK=1 for equal spacing) 
total number of normal grid points 


DHPTOL 


I ACC 


ITMAX 


IW 


Line 11 

Values 

DFPTOL 

DHPTOL 

I ACC 

ITMAX 

IW 


corresponding to line headings above 
convergence criterion for delta-F -prime 
convergence criterion for delta-H-prime 
order of accuracy, 2 or 4 
maximum number of iterations 

0 to neglect transverse curvature 

1 to include transverse curvature 


Line 12 
DFE 


DHE 


IADD 


VELEDG 


NX1 


Line 13 

Values corresponding to line headings above 

DFE criterion on F at BL edge for adding points (if IADD«1) 

DHE criterion on H at BL edge for adding points (if IADD=1) 

IADD 0 — > do not check to see if addition of normal grid points is required 

1 — > add normal grid points if required (test based on DFE, DHE) 

VELEDG value of F to be used in defining edge of BL 

NXl no. of steps (= NX-1, where NX is the total number of streamwise points) 


Line 14 
IYINT 


KODAMP 


KODPRT 


KODVIS 


KTCOD 


Line 15 

Values corresponding to line headings above (for transition/turbulence model) 
normal intermittency function set to 1 
normal intermittency function from an equation 
local values used in equation for damping 
wall values used in equation for damping 
constant PRT 

Rotta distribution for PRT 
tabular input for PRT 
mixing length model 
two-layer eddy viscosity model 
transition extent from equation 
transition extent from specified TLNGTH 


IYINT 

1 — > 
2 — > 

KODAMP 

1 — > 
2 — > 

KODPRT 

1 — > 
2 — > 
3 — > 

KODVIS 

1 — > 
2 — > 

KTKOD 

1 — > 
2 — > 
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Line 16 
SMXTR 


SST 


TLNGTH 


PRT 


NXLIM 


Line 17 

Values corresponding to line headings above {for transition/turbulence model) 
SMXTR critical vorticity Reynolds number for transition onset estimation 
SST x-location of transition (either one of SMXTR or SST will take effect) 

TLNGTH ratio of x at transition zone end to x at transition zone beginning 
PRT turb. Prandtl number 

NXLIM — > stop march at I=NXLIM (NXLIM<or=NX) 


Line 18 
RSTAR 


I0RD1 


I0RD2 


ITEMAX PRL 


Line 19 

Values corresponding to line headings above 
RSTAR gas constant, dimensional 

for air RSTAR=1716 ft*2/sec*2 degR or 286.96 m A 2/sec*2 K 
I0RD1 locations i.LE.iordl will have first order streamwise gradients 

I0RD2 locations i.GE.iord2 will have second order streamwise gradients; 

the intermediate locations will have first-order/second-order blend 
program checks for I0RD1 > 1; IORD2 > I0RD1 
ITEMAX maximum number variable entropy iterations 
PRL Prandtl number, laminar value 


Line 20 

STEP SIZES (NX1 VALUES) 


Line 21 

NXl values of step sizes, (SS (I) , 1=1, NX1) ,- values separated by a space or comma 
Line aa 

WALL OR PROFILE PRINT FLAGS (0,1 or 2), NXl VALUES 
Lines aa+1, aa+2, ... 

Print flags corresponding to each step, NXl values separated by space or comma 

0 — > no wall or profile print 

1 — > wall print 

2 — > wall and profile print 


Line bb 

OUTPUT FOR PLOTS, NUMBER OF VARIABLES FOLLOWED BY VARIABLE NUMBER CODES 
Lines bb+1, bb+2, ... 

number codes for variables to be output, see list below of variables 
that can be output; maximum 10 number codes only; first value is the 
number (IVL) of variables to be output, followed by IVL number of 
variable number codes separated by a space or comma 
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Line cc 
NUMBER 


L {INVISCID INPUTS) 


Line cc+1 

Values corresponding to line headings above 
NUMBER no. of points of inviscid data 

L order of interpolation, 1 for linear, 2 for quadratic, etc 


Line cc+2 

XE RADE SE 


PESE TWSE QESE WWSE 


Lines cc+3, cc+4, ... 

values corresponding to line headings above 


(NUMBER number of lines) 


XE axial length, m (ft) 

RADE body radius, m (ft) 

SE body surface arc length, m (ft) 

PESE BL edge pressure, Pa (lb/ft2) ( 

TWSE wall temperature, deg K (deg R) , used if IWALL-1 

QESE wall heat flux, W/m2 (Btu/ft2-s), used if IWALL-2 

WWSE wall mass flux, Pa-s/m (lb-s/ft3) 


Line dd *** lines below required only if K0DPRT=3 *** 
NUMB1 


Line dd+1 

number of values to be read in for GLAR, PRTAR table 


(NASA TM-83207 for details) 


Line dd+2 

GLAR PRTAR 


Lines dd+3, dd+4, . . . , 

NUMBl values of (GLAR, PRTAR) pairs, one pair to a line 
*** lines above required only if KODPRT=3 *** 


Line ee *** lines below required only if IENTRO-2 *** 

NS (NO. OF CURVED SHOCK COORDINATES FOR VARIABLE ENTROPY 


CALCULATION) 


Value corresponding to line heading above (number of curved shock co-ordinates) 
Line ee+2 

RRS ZZS DRSDZS 

Lines ee+3, ee+4, ... 

NS values of groups of (RRS, ZZS, DRSDZS) , one group to a line 
describing the shape of the shock line defined as below: . 

RRS radial location of a point on the shock curve, dimensions 

Z2S axial location of a point on the shock curve, dimensional 

DRSDZS derivative of RRS with respect to ZZS 

*** lines above required only if IENTR0=2 *** 
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APPENDIX C. INPUT PARAMETERS FOR BL2D 


AK 

AMACH 

CONVE 

DFE 

DFPTOL 

DHE 

DHPTOL 

DRDZS 

GAM 

GLAR 

IACC 

I ADD 

IBODY 

IENTRO 

IFS 

IFT 

I GAS 

IORD1 

IORD2 

ITMAX 

ITEMAX 

IUNIT 

IVL 

IW 

I WALL 
IYINT 
J 

KODAMP 

KODPRT 

KODVIS 

KTCOD 

L 

NS 


stretching parameter k s for inner distribution of normal grid; k s > 1 
mach number, M 0 0 

convergence criterion on u € for variable entropy iteration 

criterion on F' at boundary-layer edge for adding points (used only if 

IADD=l) 

convergence criterion on 6F' 

criterion on H' at boundary- layer edge for adding points (used only if 
IADD=l) 

convergence criterion on 8H' 

derivative of shock curve radius with respect to axial length 
ratio of specific heats, 7 

array {z/6) that corresponds to PRTAR array (see NASA TM-83207 for 
details) 

order of accuracy, 2 or 4 

0, no addition of normal grid points; 1, add if required 

1 for flow with a stagnation point; 2 for flow without stagnation point 
1 for constant entropy; 2 for variable entropy 
1 for input of PTS, TTS; -1 for input of PFS, TFS 

0 for locally similar solution; 1 for non-similar solution (recommended) 

1 for Sutherland law for viscosity; 2 for power law; if IGAS is 2, check 
power law coefficients vis2cl , vis2c2 in s/r ref 

value such that when I .LE . IORD1, streamwise gradients are computed 
to first order; I0RD1 > 1 

value such that when I . GE . IORD2, streamwise gradients are computed 
to second order; IORD2 > IORD1 
maximum number of iterations 

maximum number of iterations in the variable entropy loop 
0 for British units; 1 for SI units 

number of variables to be output from list in appendix D 
0 to neglect transverse curvature; 1 to include transverse curvature; select 
0 for internal flows 

wall bounday condition: 0 for adiabatic; 1 for temperature specified; 2 for 
heat flux specified 

1, normal intermittency function set to 1; 2, function from equation (see 
NASA TM-83207 for details) 

0 for two-dimensional; 1 for axisymmetric 

1, local values used in equation for damping; 2, wall values used 

1, constant PRT; 2, Rotta distribution for PRT; 3, table input for PRT 

1, mixing length model; 2, two-layer eddy viscosity model 

1, transition extent from equation; 2, from specified TLNGTH 

order of interpolation: 1 for linear; 2 for quadratic, etc. 

number of points that define the shock curve (used only if IENTRO = 2) 
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NUMBER 

NUMB1 

NX1 

NXLIM 

NZ 

NZI 

PESE 

PFS 

PHI I 

PRL 

PRT 

PRTAR 

PTS 

QESE 

RADE 

RRS 

RSTAR 

SE 

SMXTR 

SS 

SST 

TFS 

TLNGTH 

TTS 

TWSE 

VELEDG 

WAVE 

WWSE 

XE 

ZI 

ZMAX 

ZZS 


number of input points that define the boundary-layer edge conditions 
number of values of (GLAR, PRTAR) pairs 

number of steps ( = NX - 1, where NX is the total number of streamwise 
points) 

stop march at I = NXLIM (NXLIM. LE. NX ) 

total number of normal grid points, ke 

number of mesh points in the inner distribution, ki 

boundary-layer edge pressure, P e (N/m 2 or lb/ft 2 ) 
free-stream static pressure, Poo (N/m 2 or lb/ft 2 ) 
opening angle of body at s* = 0, deg 
Prandtl number, Np r , laminar value 
Prandtl number Np r ,t, turbulent value 

array of PRT values input in tabular form (for KODPRT = 3) 

free-stream total pressure, Po (N/m 2 or lb/ft 2 ) 

wall heat flux q* w \ used if I WALL * 2 (W/m 2 or Btu/ft 2 -sec) 

body radius (m or ft) 

shock curve radius (m or ft) 

gas constant (m 2 /(sec 2 K) or ft 2 /(sec 2o R)); for air RSTAR = 286.96 or 
1716.0 

body surface arc length; s* at input data points (m or ft), 
critical vorticity Reynolds number for transition onset estimation 
array of length NX1 containing step size values 

s* location of transition (to be set as a large value if SMXTR criterion is 
to be used) 

free-stream static temperature, (K or °R) 

ratio of s* at end of transition zone to s* at beginning of transition zone 
(for KTKOD = 2 ) 

free-stream stagnation temperature, To (K or °R) 
wall temperature, used if IWALL = 1 (K or °R) 
value of F to be used to define edge of boundary layer 
shock- wave angle at s* = 0 (input 0 if no shock) 
wall mass flux w w * (Pa-s/m or lb-s/ft 3 ) 
axial length (m or ft) 

normal coordinate value at the end of inner distribution, (, 
the maximum value of the transformed normal coordinate, 
shock curve axial location (m or ft) 
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APPENDIX D. NUMBER CODES FOR 
OUTPUT OF VARIABLES IN BL2D 


Number 

BL2D 

Description of the output variable 

code 

variable 

(output to file fort . 7) 

1 

betah 

A * 

pressure gradient parameter, 2£^- 

2 

cfe 

skin-friction coefficient based on edge conditions, Cy e 

3 

cfw 

skin-friction coefficient based on wall conditions, Cf <w 

4 

disp 

displacement thickness, 6 * (m or ft) 

5 

dpeds 

pressure gradient, d/ds* of P e 

6 

dsmxo 

d/ds* of maximum vorticity Reynolds number 

7 

dteds 

temperature gradient, d/ds* of T e 

8 

dueds 

velocity gradient, d/ds* of u* 

9 

error 

convergence parameter, dF'/F' at the wall 

10 

form 

ratio of displacement thickness to momentum thickness 

11 

hd 

heat transfer coefficient, h 

12 

itro 

number of variable entropy iterations 

13 

ame 

edge Mach number 

14 

amues 

edge viscosity, /i*, dimensional 

15 

it 

number of iterations 

16 

anste 

Stanton number based on edge condition 

17 

anstw 

Stanton number based on wall condition 

18 

anue 

Nusselt number based on edge condition 

19 

anuw 

Nusselt number based on wall condition 

20 

eps 

1 / y/ ■ -R^ref 

21 

p20 

total pressure at boundary-layer edge, nondimensional 

22 

pes (i) 

boundary-layer edge pressure at location i (N/m 2 or lb /ft 

23 

qsd 

heat transfer at the wall (W/m 2 or Btu/ft 2 -sec) 

24 

roes 

edge density (kg-sec 2 /m 4 or lb-sec 2 /ft 4 ) 

25 

redelt 

Reynolds number based on local displacement thickness 

26 

res 

local Reynolds number 
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Number 

BL2D 

Description of the output variable 

code 

variable 

(output to file fort . 7) 

27 

rethet 

Reynolds number based on local momentum thickness 

28 

rftrue 

recovery factor 

29 

rad (i) 

body radius (m or ft) 

30 

smxp 

maximum vorticity Reynolds number 

31 

rshk 

local radius of shock wave (m or ft) 

32 

rvwal 

suction or blowing mass flux normalized with edge value 

33 

rvwald 

dimensional mass flux at the wall (Pa-s/m or lb-sec/ft ) 

34 

x (i) 

boundary-layer surface coordinate s* (m or ft) 

35 

swang 

local shock-wave angle in degrees 

36 

taud 

wall shear stress (N/m 2 or lb/ft 2 ) 

37 

tes 

dimensional edge temperature, T e 

38 

theta 

momentum thickness (m or ft) 

39 

t rf ct 

intermittency distribution 

40 

twbttl 

T w /To 

41 

ues 

dimensional edge velocity, u* (m/sec or ft/sec) 

42 

utau 

friction velocity, yfrfp at wall (m/sec or ft/sec) 

43 

vw 

transformed wall normal velocity at wall, w w 

44 

alphah 

(7 - 1)M* 

45 

xi (i) 

{ 

46 

ye 

boundary-layer thickness (based on VELEDG), t* (m or ft) 

47 

jpoint 

normal grid index k, where turbulence model inner law = outer 
law 

48 

xa (i) 

axial coordinate of body, x* (m or ft) 

49 

zshk 

axial coordinate of shock wave (m or ft) 
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APPENDIX E. DESCRIPTION OF MAIN 
COMPUTER VARIABLES IN BL2D 


BL2D variable 

adl, ad2, ad3 
asl, as2, as3 
ak 

al, alp 
all, allp 
al2, al2p 
amach, amache 
amu, amue 
amues 
amuf s 
amuref s 
amurs 

anste, anstw 
anue , anuw 
bet, bet 2 
betah 

cfe, cfw 
conve 

delfp, delhp 

dfptol, dhptol 

disp 

dpeds 

drdzs 

ds 

dsmxo 
dxi, dxil 


Description 


01 , 02 , 03 ; coefficients used in streamwise differencing d/d£ 
coefficients used in streamwise differencing d/ds* 
normal grid stretching parameter, k s 

IJ 

h,l[ 

hA 

Moo, M e 
M, Me 

edge viscosity, /i* , dimensional 
Aoo 

n* at T re f 

Sutherland viscosity law constant 
Stanton number based on edge or wall condition 
Nusselt number based on edge or wall condition 
7/(7- 1), (7- l)/2 
pressure gradient parameter, 

skin-friction coefficient based on edge or wall condition 
convergence criterion for variable entropy iteration 
8F',8H' 

convergence criteria on 8F', 8H' 
displacement thickness, 8* (m or ft) 
pressure gradient, 
slope of the shock curve 
step size in s* 

d/ds* of maximum vorticity Reynolds number 
step sizes in & - &_ lf &_i - f,_ 2 
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BL2D variable 


Description 


dsi, dsil 

step sizes in s* 

eps 

1 / \/ R e ie{ 

f,fl,f2 

F at i, i — 1, i — 2 

form 

ratio of displacement thickness to momentum thickness 

fp, fpl, fp2 

F' at i, i — 1, i — 2 

g, gam 

7 

glar 

array of z* / 6* values for which Np Ttt values are input 

h,hhl, hh2 

H at i, i — 1, i — 2 

hd 

heat transfer coefficient, h 

hp, hpl , hp2 

H' at i, i — 1, i — 2 

i, il 

i, i — 1 

iacc 

order of accuracy 

iadd 

flag for addition of normal grid points 

ibody 

flag for type of flow at i = 1 

ientro 

flag for variable entropy calculation 

ift 

flag for nonsimilar solution 

if s 

flag for input of free-stream conditions 

igas 

flag for type of viscosity law 

iordl, iord2 

i indices to specify order of accuracy of the streamwise gradient 

it, itmax 

number of iterations and maximum number of iterations 

ite, itemax 

number of variable entropy iterations and its upper limit 

itrans 

flag to denote laminar, transitional, or turbulent regime (0, 1, or 2) 

itro 

number of variable entropy iterations 

iunit 

flag for choice of units, US or SI 

iw 

flag for transverse curvature term 

iwall 

flag for energy boundary condition at the wall 

iyint 

normal intermittency function 

j 

flag for two-dimensional or axisymmetric flow 

k, kl 

normal grid index k, k — 1 

nmax 

greater value of (nzm, nxm) 

ns, nsm 

number of points defining the shock curve, upper limit for ns 
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BL2D variable Description 


num 

number of input inviscid data points 

numb 1, numb lm 

number of input Np T i values, upper limit for numbl 

nx, nxl 

number of streamwise grid points, number of steps 

nxlim 

streamwise grid index at which to stop march, nxlim < nx 

nxm 

maximum array dimension for streamwise grid 

nz, nzm 

number of normal grid points, upper limit for nz 

nzi 

number of normal grid points in the inner part of a two-part 
distribution 

oz 

\/2? 

plO 

total pressure at boundary-layer edge for ientro = 1, 
non-dimensional 

p20 

total pressure at boundary- layer edge for ientro = 2, 
non-dimensional 

pe 

static pressure at boundary-layer edge at grid locations 

pes 

static pressure at boundary-layer edge, dimensional 

pese 

static pressure at boundary-layer edge at input data points, 
dimensional 

pfs 

free-stream static pressure, dimensional 

phii 

opening angle of body at s* = 0, deg 

prl,prt 

Prandtl number, laminar or turbulent 

prtar 

input array of turbulent Prandtl number values 

pts 

free-stream total pressure, dimensional 

qsd, qws 

heat transfer at the wall (W/m 2 or Btu/ft 2 sec) 

qwse 

heat transfer at the wall at input data points (W /m 2 or 
Btu/ft 2 -sec) 

rad 

body radius, dimensional 

rade 

body radius at input data points 

rede It 

Reynolds number based on local displacement thickness 

refs 

free-stream Reynolds number 

reref 

Re re f 

rethet 

Reynolds number based on local momentum thickness 

roe 

edge density, nondimensional 
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BL2D variable 


Description 


rof s 

free-stream density, dimensional 

rrs 

input shock curve radius, dimensional 

rstar 

gas constant, dimensional 

ss 

array of step sizes 

sst 

transition location 

taud 

wall shear stress (N/m 2 or lb/ft 2 ) 

tio 

boundary-layer edge total temperature, nondimensional 

te 

boundary-layer edge temperature, nondimensional 

tf 

free-stream temperature, nondimensional 

tf s 

free-stream temperature, dimensional 

theta 

momentum thickness, dimensional 

tlngth 

ratio of s* at the end of transition zone to s* at the beginning of 
transition zone 

tts 

total temperature, dimensional 

tws 

wall temperature, dimensional 

twse 

wall temperature at input data points, dimensional 

ue 

edge velocity, nondimensional 

uee 

edge velocity for ientro = 2, nondimensional 

uf s 

free-stream velocity, dimensional 

veledg 

value of F used to define edge of boundary layer 

vis2cl , vis2c2 

viscosity coefficients for igas = 2 

w 

transformed normal velocity 

wave 

shock wave angle at s * — 0 

wws 

wall-normal velocity, dimensional 

wwse 

wall-normal velocity at input data points, dimensional 

x (i) 

body surface coordinate s*, dimensional 

xa (i) 

axial length, dimensional 

xae (i) 

axial length at input data points, dimensional 

xe (i) 

body surface coordinate values at input data points, dimensional 

xi 

t 

z 

transformed normal grid array, ( 
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BL2D variable 


Description 


zi 

zmax 

zshk 


value of ( that corresponds to the limit of the inner distribution 

maximum value of ( 

axial coordinate of shock wave (m or ft) 
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